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ABSTRACT 

Several localized versions of the ensemble Kalman filter have been proposed. Although tests 
applying such schemes have proven them to be extremely promising, a full basic understanding 
of the rationale and limitations of localization is currently lacking. It is one of the goals of this 
paper to contribute toward addressing this issue. The second goal is to elucidate the role played 
by chaotic wave dynamics in the propagation of information and the resulting impact on forecasts. 
To accomplish these goals, the principal tool used here will be analysis and interpretation of 
numerical experiments on a toy atmospheric model introduced by Lorenz in 2005. Propagation 
of the wave packets of this model is shown. It is found that, when an ensemble Kalman filter 
scheme is employed, the spatial correlation function obtained at each forecast cycle by averaging 
over the background ensemble members is short ranged, and this is in strong contrast to the much 
longer range correlation function obtained by averaging over states from free evolution of the 
model. Propagation of the effects of observations made in one region on forecasts in other regions is 
studied. The error covariance matrices from the analyses with localization and without localization 
are compared. From this study, major characteristics of the localization process and information 
propagation arc extracted and summarized. 



1. Introduction 

In an ensemble Kalman filter (EnKF), at any given 
analysis time t = nT, the estimated state of the atmo- 
sphere and the corresponding error covariance reflecting 
uncertainty in the state estimation are represented by a fi- 
nite ensemble of global system states (e.g., Evensen 1994; 
Burgers ct al. 1998; Houtckamer and MitchcU 1998). Each 
ensemble member is then integrated forward in time using 
a physical model of the atmosphere, thus creating an en- 
semble of forecasts at the next analysis time, t = {n+ l)T. 
By suitable incorporation (assimilation) of new measure- 
ments, the data assimilation process creates a new ensem- 
ble of system states that reflects the most probable atmo- 
spheric state and its uncertainty based on the available 
combined knowledge contained in the forecast ensemble at 
time {n + 1)T and the new measurement data. The process 
then repeats at the t = {n + 2)T cycle, and so on. 

Several localized versions of EnKF have been proposed 
(e.g., Houtckamer and MitchcU 2001; Hamill et al. 2001; 
Anderson 2001, 2007; Ott et al. 2004; Hunt et al. 2007). 
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In this paper, we concentrate on one localization scheme 
called the local ensemble Kalman filter (LEKF; Ott ct al. 
2004) and a subsequent computationally more efficient ver- 
sion, the local ensemble transform Kalman filter (LETKF; 
Hunt et al. 2007). We expect that similar results would 
be obtained using any other of the available EnKF local- 
ization schemes. In LEKF and LETKF, local regions sur- 
rounding each gridpoint location are defined, the analysis 
is performed in each local region, and the results are com- 
bined to form a global analysis ensemble. The motivation 
for using localization is that it is expected to work well 
with drastically fewer ensemble members than would be 
the case if localization were not employed. In fact, without 
localization, the computational requirements for a mean- 
ingful EnKF analysis arc so vast that they make the ap- 
proach completely infeasible in typical weather forecasting 
settings. On the other hand, the use of localization has 
been shown to be successful and eminently feasible for both 
toy models and real atmospheric models (e.g., Houtckamer 
and Mitcheh 2005; Houtekamer et al. 2005; Whitaker et al. 
2008; Szunyogh et al. 2008). 

One goal of this paper is to study the basic proper- 



ties and rationale for localization. To do this, we will 
employ numerical experiments on a simple model system 
of equations introduced by Lorenz (Lorenz 2005). The 
model is simple and small enough that (unlike a real op- 
erational model) we can employ an EnKF without local- 
ization. We will compare the results obtained with the 
localized LETKF scheme with those obtained without lo- 
calization [henceforth referred to as ensemble transform 
Kalman filter (ETKF); ETKF was introduced in Bishop 
et al. (2001)]. By doing this, we will be able to study some 
aspects relevant to the issues of localization, such as what 
determines a good size of the localization region, what dy- 
namics justifies (or not) localization, and so on. 

The second goal of this paper is to study the propa- 
gation of information in a spatiotemporally chaotic sys- 
tem with wave-like dynamics that are similar to that of 
weather. In particular, we use the traditional approach 
of discussing the propagation of information based on the 
study of phase and group velocities of the waves in the 
model (e.g., Charney 1949; Persson 2000; Szunyogh et al. 
2002; Zimin et al. 2003). An alternative approach that has 
recently been gaining increasing attention is based on us- 
ing tools of probability theory — most important, measures 
based on relative entropy of information (e.g., Kleeman 
2007). The relation between localization and the advec- 
tion of information from the observation location by the 
model dynamics has been studied using the hierarchical 
ensemble filter by Anderson (2007). 

The organization of this paper is as follows: In section 
2, wc introduce the model system (Lorenz 2005) that we 
employ, and we investigate some of its salient properties. 
In section 3, we investigate correlations of the model. In 
section 4, we explain the forecast and localization proce- 
dures used in the ETKF and LETKF schemes. In section 
5, we study correlations of the background ensemble of an 
ETKF as a function of spatial distance. In section 6, we 
compare covariance matrices obtained using an ETKF and 
an LETKF. In section 7, we show how observations taken 
in a small region at a certain time affect the forecast at 
other spatial points in the future. Section 8 provides fur- 
ther discussions and summarizes our main conclusions. 

2. Model 

Lorenz (2005) discusses three one-dimensional toy mod- 
els that incorporate many features shown in real atmo- 
spheric dynamics and in global numerical weather predic- 
tion models. The first model (Lorenz model 1) was origi- 
nally introduced in Lorenz (1996) and Lorenz and Emanuel 
(1998). This model has become the standard model of 
choice for the initial testing of EnKF schemes. The popu- 
larity of the model is in part due to the similarity between 
the propagation of uncertainties (forecast errors) in Lorenz 
model 1 and global models in the midlatitudc storm-track 



regions. In particular, the errors are propagated by disper- 
sive waves whose behavior is similar to that of synoptic- 
scale Rossby waves, and the magnitude of the errors has 
a doubling time of about 1.5 days (where the dimension- 
less model time has been converted to dimensional time by 
assuming that the characteristic dissipation time scale in 
the real atmosphere is 5 days; see Lorenz 1996). Lorenz 
model 2 adds the feature of a smooth spatial variation of 
the model variables that resembles the smooth variation 
of the geopotential height (streamfunction) at the synoptic 
and large scales in the atmosphere. Lorenz model 3, the 
most refined and "realistic" of the three models in Lorenz 
(2005), adds a rapidly varying small-amplitude component 
to the smooth large-scale flow, mimicking the effects of 
small-scale atmospheric processes. In our following simu- 
lations, we use Lorenz model 3. 

The equation of evolution of the scalar state variable 
Zn at position n is the following: 

dZr^ldt = + 62[r,r]i,„ + c[y,x]i,„ 

- x„ - 6y„ + (1) 

where n is a spatial variable (n = l,2,...,Af), periodic 
boundary conditions are employed (Zjv+i = ^i), and the 
X and Y vectors arc defined as 

= 51 - (2) 

Yn ~ Zn — Xn, (3) 

with a = (3/2 + 3)/(2/3 + 4/), 
/3 = (2/2 + l)/(/4 + 2/2). 

The prime notation on S' signifies that the first and the 
last terms in the summation are divided by 2. The bracket 
of any two vectors X and Y is defined as 

J ,/ 

[^,Y]K,n= ^ ' ^ '{~Xn^2K-iYn-K~j 

+ X„_K+,_,r„+K+,)//^2 (4) 

when K is even, and S' is replaced by S when K is odd; 
J = K/2 when K is even, and J = (K — l)/2 when K 
is odd. The parameter values used throughout this paper 
are N = 960, /C = 32, & = 10, c = 2.5, F = 15, and / = 12. 
We solve Eqs. (l)-(4) using a Runge-Kutta scheme with a 
time step sufficiently small to resolve the fast time scale of 
the model. We find that these equations lead to a state 
profile that shows small-scale activity K„ superposed on a 
large-scale smooth component X„, thus crudely mimicking 
the multiscale dynamics of a real atmospheric system. 

Figure 1 obtained from a numerical solution of Eq. (1) 
shows how states evolve in time. The horizontal axis is 
the spatial location, and the vertical axis is the time. We 
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Fig. 1. Time evolution of state Zn of Lorenz model 3. 
In this figure, 6 h are defined to be 0.05 time units in 
the model equations considering the damping time of the 
model. 



Fig. 2. Time evolution of the wave packet envelopes ex- 
tracted by the envelope extraction method introduced in 
Zimin et al. (2003). Wave numbers of only 7 and 8 were 
used to smooth out the envelopes. 



define 6 h to be 0.05 time units in the time evolution equa- 
tion following Lorenz [Lorenz justifies this correspondence 
between t in Eq. (1) and pseudohours by considering the 
damping time of the model] . The colors in Fig. 1 represent 
the values of the state variables Zn- Red corresponds to 
high values, and blue corresponds to low values. The figure 
shows that waves of approximately 7 spatial periods prop- 
agate to the left. This observed structure mimics what is 
seen at large scale in the atmosphere where Rossby waves 
play a prominent role. 

The structure found in Fig. 1 also suggests that an in- 
sight can be gained by representing the field as a modu- 
lated sinusoidal wave. For this purpose, modulating en- 
velopes were extracted from the states at each time using 
the method introduced in the paper by Zimin et al. (2003). 
The method of Zimin et al. for extracting modulating en- 
velopes is explained briefly in appendix A. Based on the 
observed structure in Fig. 1, we computed envelopes using 
wavenumbers of only 7 and 8 to show smooth flow of the 
envelope (see appendix A). Figure 2 depicts the envelope 
amplitude in time and space for the same numerical run as 
in Fig. 1. From Fig. 2, we see that wave packets propagate 
to the right, in contrast to the leftward propagation of the 
individual troughs and ridges seen in Fig. 1. That is, in 
some appropriate sense, the wave turbulent state consists 
of waves that have phase velocities that move to the left 
(as in Fig. 1) and group velocities that move to the right 
(as in Fig. 2). This situation is similar to that of Rossby 
waves in the atmosphere, whose phase velocity is always 
westward relative to the mean zonal flow and whose group 
velocity can be either eastward or westward. 



3. Correlation structure of the model 

We calculated the correlations Cm (j) between values at 
two different locations separated by m grid points in space 
and by t in time as follows: 

r, , , _ ( [^n(t) - (^n(^)) ] [^n+n,(t + t) - (Z„(t)) ] ) 

([Z„(t)-(Z„(t))]^) 

(5) 

where angle brackets denote an average over n and t. The 
resulting plots, Cm{T) (Fig. 3) and Cm(0) (Fig. 4), show 
that wave dynamics plays a central role in the spatiotem- 
poral evolution of the correlation, as demonstrated through 
three properties: (1) the spatial correlation has the struc- 
ture of a wave packet with carrier wavcnumber 7 and a 
packet envelope that decreases monotonically with increas- 
ing |m| (Fig. 4), (2) the wave structure shifts to the left 
with a phase velocity of about -0.77 grid points per hour 
(Fig. 3), and (3) the packet envelope shifts to the right at 
a rate that we identify as a group velocity (e.g., see the 
bright yellow spots along the red troughs in Fig. 3). Under 
the assumption of ergodicity, the correlation in Eq. (5) can 
be regarded as describing an ensemble drawn from the dis- 
tribution that defines the climate of the model. In section 
5, we will consider another type of correlation function that 
can be regarded as describing an ensemble drawn from a 
distribution of short-term forecast uncertainties associated 
with a fairly dense observational network. 

To make property (3) more transparent, we extracted 
the envelopes at each temporal distance t from Fig. 3 us- 
ing the envelope extraction method in Zimin et al. (2003) 
with wavenumbers of 6, 7, and 8. We obtained the results 
depicted in Fig. 5. It is seen from Fig. 5 that, similar to 
Fig. 2, the envelope in m-r space propagates to the right 
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Fig. 3. Correlations Cm{T) between state variables at two 
different points, {n,t) and (n + m,t + r) in Fig. 1. The 
values were averaged over all of the space-time locations 
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Fig. 4. Correlations C„i(0); i.e., the correlation values of 
Fig. 3 with time difference r set to 0. 



at a group velocity of about 1.37 grid points per hour. 

4. Forecasting and localization 

In our numerical assimilation experiments, we will em- 
ploy "perfect model" simulations of forecast/analysis cy- 
cles. That is, we carry out a time evolution run of our 
model, Eqs. (l)-(4), which we regard as simulating a "true" 
atmospheric state evolution that we wish to analyze and 
forecast. We then take simulated measurements of this 
true state Zn by adding uncorrelated noise of Gaussian dis- 
tribution to the true state variable Z„ at each grid point 
and at each analysis time. The mean of the noise is 
and the standard deviation is 0.3, which is about the size 



Fig. 5. Envelopes extracted from correlations in Fig. 3 at 
each time r using wave numbers only from 6 to 8 to smooth 
out the result. 



of the small-scale activity Y in the model. The standard 
deviation of the noise can be also compared with the clima- 
tological standard deviation of the state variable Zm which 
is 4.67 (a noise standard deviation of 0.3 is used through- 
out this paper). Using these simulated measurements, we 
do analyses employing the same model, Eqs. (l)-(4), as we 
use to generate the true state evolutions [i.e., Eqs. (l)-(4)]. 

In our ETKF implementation, we use the same tech- 
nique as described in Hunt et al. (2007), but without em- 
ploying localization. To faithfully represent the full system 
state and its covariance in the absence of localization, the 
number of ensemble members k that we use in the ETKF 
must be at least on the order of the number of growing 
local Lyapunov exponents, which increases as the number 
of grid points of the model increases in general. Thus, we 
need many ensemble members when there are many grid 
points in the model. We integrate each ensemble mem- 
ber during the evolution phase of the forecast cycle, and 
the required computing power to do this scales like k. In 
carrying out our ETKF procedure, it is also necessary to 
invert a kxk matrix, which requires a number of computer 
operations on the order of k^. Hence, if we can do compu- 
tations in local regions with much fewer grid points than 
in the global region, and thus with correspondingly much 
fewer ensemble members, then the required computational 
power is greatly reduced. For a more detailed discussion of 
the computational cost of LETKF, see Hunt et al. (2007) 
and Szunyogh et al. (2008). 

If the errors in the state estimates at two different grid 
points that are far from each other are independent, then 
we might be able to compute analysis ensembles for each 
location while ignoring the other. This qualitative type of 
consideration is what motivates localization of the analy- 
sis. A brief description of the localization process used in 
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LETKF and implemented on our one-dimensional model, 
Eq. (1), is the following: First, choose an appropriate spa- 
tial size L for the localization. At each grid point z, a local 
patch is assigned as consisting of the 2L + I grid points 
{i - L,i - L + 1,. . . ,i + L - 1,1 + L) for i = 1, 2, . . . , iV. 
We then compute the analysis ensemble in each local patch 
using observations in the local patch. Then, we construct 
the global analysis ensemble that combines all of the local 
analysis ensembles by taking weighted averages of analy- 
sis values from each local patch. Ott et al. (2004), Hunt 
et al. (2007), and appendix B explain how to do the steps 
in this procedure in detail. How accurate is this localiza- 
tion procedure? The supposed main requirement is that 
observations outside the local patch should not affect the 
information from the local patch that is used in the evolu- 
tion of the ensemble to the next analysis cycle. We discuss 
this issue for Lorcnz model 3 in what follows. Note that, 
because we investigate a univariate model, the issue of bal- 
ance does not arise. The issue of balance and its interaction 
with localization has been discussed by Kepcrt (2009). 

5. Spatial correlations of the background ensemble 

As we have pointed out. Fig. 4 shows a long-range cor- 
relation essentially extending over the entire length of the 
simulation region. Thus, one might think that we will lose 
a lot of information that comes through correlations from 
observations outside the local patch when we localize the 
analysis scheme. As we shall subsequently argue, this is 
not the case. In particular, we emphasize that it is the 
covariance matrix of the background ensemble around the 
ensemble mean that is used in the analysis, not the covari- 
ance matrix of the climatological distribution of true states 
themselves. 

To formulate a correlation relevant to what is used in 
our analyses, we first carry out a perfect-model simula- 
tion of a forecasting situation. In particular, we begin by 
preparing a background ensemble Zn^\t),k = 1,2,..., 960, 
by running 1000 forecast cycles (we use k and t to denote 
an ensemble member and cycle time, respectively). Dur- 
ing each cycle, we saved the background ensemble (t) 
and took simulated observations with noise of uncorrelated 
Gaussian distribution at each grid point at the analysis 
time, and then we updated the ensemble with ETKF using 
a multiplicative covariance inflation factor of 1.13 [see step 
5 of appendix B; the factor 1.13 gives the minimum rms 
error of the resulting state estimate — see Anderson and An- 
derson (1999) for a general explanation of covariance infla- 
tion] , and we then evolved the true state and the ensemble 
from ttot + Qh (except where otherwise stated, a 6-h cycle 
time and 1.13 multiplicative covariance inflation factor are 
used throughout this paper). Using these parameters, we 
calculated the following ensemble-based spatial correlation 
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Fig. 6. Correlation of the background ensemble of 
ETKF with 6-h cycle time vs spatial distance m. The 
values were averaged over 1000 forecast cycles. 

in the background: 

with A = ( [Z,l^-)(t) - Z„(t) ] - Zn+m{t) ] 

i? = ([zW(t)-Z„(t)]2),, 

where (•)fc denotes an average over ensemble members fc, 
{■)n,t denotes an average over space n and time t, and 
Zn ~ { zli'^ )k- Figure 6 shows the result. As can be seen, 
the magnitude of the correlation is less than 3% of its peak 
value for \m\ > 46, and the half-width of the strong peak 
around zero distance is roughly on the order of 10 for 6-h 
cycle time simulation. Thus, the correlations of deviations 
of ensemble members from the ensemble mean will be neg- 
ligible when 2 grid points are separated by more than 46 
grid points. 

To analyze further the result shown in Fig. 6, we note 
that the background correlation is determined by the struc- 
ture (correlation) of the analysis error and the model dy- 
namics that spreads the effects of the analysis error in space 
and time. The shape of the analysis correlation function 
is very similar to that of the background correlation func- 
tion and is only weakly dependent on the magnitude of the 
observation noise (Figs. 7a, b). A reasonable assumption 
is that the spatial expansion of the correlation pattern of 
the analysis ensemble during the model integration is con- 
trolled by the group velocity of the waves that compose 
the envelope of the analysis error pattern. Although the 
group velocity associated with these waves is not necessar- 
ily identical to the group velocity of the waves observed 
in the free run of the model, the group velocity observed 
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Fig. 7. Correlations (a) of the background ensembles 
and (b) of the analysis ensembles of ETKF with several 
different noise sizes. The values were averaged over 1000 
forecast cycles. The red, blue, brown, and black lines arc 
for standard deviations of the noise: 1.0, 0.1, 0.01, and 
0.001, respectively. 

for the free run can be considered to be an estimate of the 
maximum group velocity for the given dynamical system. 
Thus, the comparatively short range of the correlations is 
natural if one assumes that the expansion of the observed 
half-width of the correlation pattern of the analysis ensem- 
ble during the model integration cannot be larger than the 
maximum group velocity multiplied by the forecast time 
used in calculating the values Zn^\t) in Eq. (6) and iden- 
tifies the maximum group velocity (maximum information 
propagation speed) to be on the order of the speeds found 
in Figs. 1-5 (i.e., ^1 grid point per hour). 

In Fig. 8, we show results for the same quantity as 
plotted in Fig. 6, but for the cases in which the assimi- 
lation cycle times are 6 (as in Fig. 6), 24, and 48 h. The 
observations were taken only at analysis times with cor- 
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Fig. 8. Correlation of the background ensembles of 
ETKF with various cycle times vs spatial distance m. The 
values were averaged over 1000 forecast cycles. The solid, 
dash-dotted, and dotted lines are for 6-, 24-, and 48-h cycle 
times, respectively. 

responding intervals of 6, 24, and 48 h. Consistent with 
the above interpretation, we observe that the correlation 
function spreads as the cycle time is increased. 

6. ETKF and LETKF covariance matrices 

Now we compare the covariance matrices obtained from 
ETKF analyses with those obtained from LETKF analy- 
ses. We ran forecast cycles keeping two sets of ensem- 
bles, one for ETKF and one for LETKF. The number 
of ensemble members used in our ETKF implementation 
was 960, which is the total number of grid points of the 
model. Our LETKF implementation used a local patch 
size of 2 X 50 -|- 1 and number of ensemble members equal 
to 101, the same as the number of grid points in the local 
patch. The multiplicative covariance inflation factor was 
1.13 for both ETKF and LETKF, which gives the mini- 
mum rms state estimate errors for both. We computed the 
covariance matrix of the ETKF analysis ensemble and the 
covariance matrix of the LETKF analysis ensemble at sev- 
eral analysis times. Because they have nonnegligible values 
only near the diagonal, we chose the elements Pi,j of the 
covariance matrices with |j — i| < 7 and stacked them row 
by row so that elements with the same separation of indices 
j — i are aligned vertically. Figures 9a, b show the ETKF 
and LETKF covariance matrices at a representative anal- 
ysis time. The vertical axis is index i, the horizontal axis 
is j — i, and the color represents the covariance values be- 
tween grid points i and j: Pij. Red corresponds to high 
values, and blue corresponds to low values. These figures 
show similar patterns for ETKF and LETKF. Figure 10 is 
a plot of Pi versus i, where Pi is the average of Pi+j^i+j 
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Fig. 9. Covariance matrix values of analysis ensemble at 
a certain forecast time of (a) ETKF and (b) LETKF. The 
y axis is the row index i, and the x axis is the difference 
between the column indices j and i. 

over the 11 values of I for — 5 < ^ < 5. The red color 
is for ETKF, and the blue color is for LETKF. This fig- 
ure also shows remarkably similar patterns for ETKF and 
LETKF. The mean rms errors of the state estimates [de- 
fined as ( \/{{Z^~~Znj^yri) cycle, whcrc is thc mean of 
the analysis ensemble and (• • • )„ and (• • • )cycie are aver- 
ages over grid points and over forecast cycles, respectively] 
obtained from ETKF and LETKF are 0.0979 and 0.0997. 

We also studied the patch size and cycle time depen- 
dence of two global measures of the goodness of LETKF 
assimilations for 6-, 12-, 24-, and 48-h cycle times. One 
measure is the time-averaged rms error of the state esti- 
mate by averaged ensemble as defined above. The other 
measure is 
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Fig. 10. Variance values of analysis ensemble at a certain 
forecast time. The values were averaged over 11 neighbor- 
ing grid points. Red plot is for ETKF, and the blue plot is 
for LETKF. 



whcrc {■ ■ ■)i denotes an average over all of the grid points. 
We used the same number of ensemble members as the 
number of grid points in a local patch for LETKF. Figures 
11 and 12 show that the rms error and the quantity 5 be- 
come smaller as we increase the local patch size. Using the 
rms error as our figure of merit, we see that there is rela- 
tively little benefit to increasing the patch size past about 
15, although there is more substantial improvement in 5 
as the patch size is increased past 15. Thus, we see that 
the improvement in the estimate of the mean of the en- 
semble saturates more quickly with increase in patch size 
than the improvement in the estimate of the variance of 
the ensemble. From Fig. 12, we can also see that longer 
forecast cycle time makes 5 larger when the local patch size 
is smaller than about 20. However, we cannot see a clear 
pattern when the local patch size is large. Figure 13 shows 
the difference between the rms errors of LETKF and ETKF 
divided by the rms error of ETKF, a normalized rms error 
difference. 

7. Demonstration of the localized influence of ob- 
servations 

If the local patch size is large enough to encompass 
most of the correlation seen in Fig. 6, then it is to be ex- 
pected that the analysis at grid points around the middle 
of the patch does not lose too much information by ignor- 
ing observations outside the local patch. In what follows, 
we report evidence supporting this intuition. 

We now examine how observations taken at time t and 
in a small region around spatial location i affect a forecast 
at space-time point (i-l-r, i+m). To do this, we first ran 100 
forecast cycles using the ETKF and noisy observations at 
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Fig. 11. Averaged rms errors of the state estimates by 
the means of the analysis ensembles vs local patch size. 
The values were averaged over 2000 forecast cycles. The 
lines with the circles, squares, triangles, and times signs 
correspond to 6-, 12-, 24-, and 48-h cycle times. 
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Fig. 12. Values of 6 (normalized difference between vari- 
ances of analysis ensembles of ETKF and LETKF) vs local 
patch size. The values were averaged over 2000 forecast cy- 
cles. The black, red, cyan, and blue colors correspond to 
6-, 12-, 24-, and 48-h cycle times. 

all the grid points (as done in section 5), thus bringing the 
ensemble members close to the true state. Then, at current 
time t, we took observations at only five fixed contiguous 
locations, i — 2,i—l,i,i + l, and i + 2 (taking observations 
at only one location did not produce a clear result, possi- 
bly because the impact of one observation was not strong 
enough). We then used these observations to calculate an 
analysis ensemble. Then, we calculated the mean of the 
background ensemble and the mean of the analysis ensem- 
ble and evolved the two means from the current time t to 



Fig. 13. Normalized difference between rms errors of state 
estimates with ETKF and LETKF vs local patch size. The 
values were averaged over 2000 forecast cycles. The black, 
red, cyan, and blue colors correspond to 6-, 12-, 24-, and 
48-h cycle times. 

the time t + t. We denote the values of these two evolved 
means at each grid point n by Z^^ir) and Z^^{t) (b and a 
denote background and analysis, respectively). We did this 
for values of r ranging from r = h to t = 150 h in 6-h 
intervals, calculating and saving the quantities 

I?„(r) = (%„(r)-Z^+„(r))^ (8) 

We then proceed to the next time cycle t ^ t + Q h and 
repeat this process: i.e., at time t we now take observa- 
tions at all the grid points, calculate the analysis ensem- 
ble, evolve it and the true state forward to create a new 
background ensemble at the new time t ^ t + 6 h, and re- 
peat the previously described five-observation-point calcu- 
lation of Dmir). Figure 14a shows a plot of {D,n{T)) cycle, 
the average of D„i{t) over 4000 cycles. The horizontal 
axis is m, and the vertical axis is r. The color represents 
{Dmij)) cycle- Red corresponds to high values, and blue 
corresponds to low values. The white color represents val- 
ues higher than the upper limit of the color bar. The figure 
shows that observations do not affect points that are far 
from the observation points within up to about 30 h. How- 
ever, all of the points become affected by the observations 
after about 30 h. 

To investigate why all of the points eventually become 
affected after some time, we plotted the time evolution of 
the differences between the two means, 

d™(rO = Zf+„(r,)-Zf+„(r,), (9) 

at T,; = i X 6 h (i = 0, 1, . . . , 7) at a certain cycle time (see 
Fig. 15a). Figure 15a shows dmin) plotted versus m at 
6-h intervals from r = h to r = 42 h. It is seen that 
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Fig. 14. Values of (An (''')) cycie (evolution of the pertur- 
bation to the ensemble mean caused by five localized ob- 
servations). The values were averaged over 4000 cycles, (a) 
with ETKF and (b) with LETKF. The white line shows 
the group velocity of Lorcnz model 3, and the two red lines 
show the boundaries of the wedgelike region of the pertur- 
bations. 
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Fig. 15. Evolution of the perturbation to the ensemble 
mean caused by five localized observations, dm{Ti), in 6-h 
time intervals from r = h to r = 42 h, at a certain cycle 
time, (a) with ETKF and (b) with LETKF. The solid line 
shows the group velocity of Lorcnz model 3, and the two 
dotted lines show the same boundaries in Fig. 14b. 



small differences at points far from the observation points 
at r = h become amplified and dominate the dm values at 
later times. At small r, these differences are close to zero 
but are not exactly zero. We computed the rms values of 
drn{0) averaged over grid points m and cycles varying the 
number of ensemble members, where the average over m is 
taken only for m < —300 and m > 300. Figure 16 shows 
that the impact of the localized observation on points far 
from the observation point decreases as the ensemble size 
increases. So, these small impacts are likely to be caused by 
the finite size of the ensemble [see Hamill et al. (2001) and 
Anderson (2007) for detailed discussions of noisy covari- 
ance values]. As t increases, the small perturbations far 
from the observation points are chaotically amplified and 
become large. We repeated this experiment using LETKF 



in place of ETKF. The use of localization in LETKF has 
the effect of eliminating the randomlike small initial differ- 
ences c?r)i(0) at large \m\. Figures 14b and 15b show the 
results. It is seen that the effect of the localized observa- 
tions spreads to other points linearly, creating a wedgelikc 
region in our r-vcrsus-m diagram indicated by two red lines 
in Fig. 14b and two dotted lines in Fig. 15b. The propa- 
gation speed corresponding to the right boundary of the 
wedge is roughly 3 grid points per hour, and the propaga- 
tion speed corresponding to the left boundary of the wedge 
is roughly 1 grid point per hour. Thus, it appears that, 
during a 6-h cycle time, the effect of an observation will, 
on average, reach up to about 18 grid points to the right 
and up to about 6 grid points to the left. The maximum 
difference between the background mean and the analysis 
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Fig. 16. The rms values of d,„(0) divided by the size of 
the noise of observations vs ensemble size of ETKF. The 
values were averaged over 2000 cycles and over grid points 
that are more than 300 grid points away from the middle 
of the observation points. 



mean propagates at a speed that is close to the group ve- 
locity of wave packets (1.37 grid points per hour, shown by 
the white line in Fig. 14b and the solid line in Fig. 15b) 
determined from Fig. 5. 

To see the effect of observations at a space-time point 
(t, n) on a forecast at (t + T,n + m) when observations are 
taken at all of the grid points, we did an experiment similar 
to the one described above. This time, we took observa- 
tions at all of the grid points (instead of at 5 adjoining grid 
points) and calculated the value Cm(r) instead of D„i(t) 
defined in the above experiment: 



n. cycle 

(10) 

where Cnoise denotes the standard deviation of the ob- 
servation noise. Figure 17 shows C„i(t). The horizontal 
axis is m, and the vertical axis is r. The color represents 
C„j(r). Red corresponds to high values, and blue corre- 
sponds to low values. The white color represents values 
that are higher than the upper limit of the color bar. This 
plot also shows that the effect of observations propagates 
globally to the right at a speed that is close to the group 
velocity of wave packets (shown by a long white line) . The 
orientation of major axes of the blue blobs slopes to the 
left, and this slope is on the order of the phase velocity 
(-0.77 grid points per hour, shown by short white lines) 
determined from Fig. 3. 
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Fig. 17. Values of Cm (t) [correlations between perturba- 
tions at (0,n) and at (r, n -|- m) caused by observations]. 
The observations were taken at all the grid points. The val- 
ues were averaged over grid points n and forecast cycles. 
The long white line shows the group velocity of the Lorenz 
model 3, and the short white lines show phase velocity of 
the model. 



8. Conclusions 

The correlation function obtained from the time evolu- 
tion of the state of Lorenz model 3 and the envelopes ex- 
tracted from it both show that there is predominant wave 
packet propagating to the right (sections 2 and 3). We 
found that the correlation length of the deviations from the 
ensemble mean of the background ensemble of the ETKF 
at each forecast cycle with Lorenz model 3 is very much 
shorter than the correlation length of the climatological dis- 
tribution of model states (section 5). Thus, we argued that 
we do not lose much information from localization of the 
analysis if the size of the local patch is big enough to cover 
the correlation length of the deviations of the background 
ensemble. The comparison of the covariance matrices of 
the analysis ensembles with and without localization shows 
that they have similar patterns (section 6), thus providing 
strong support for the achievable accuracy of localization. 
The effect of an observation at space-time point (t, n) on a 
forecast at (t-\-T,n + m) was found to be local. In addition, 
the information obtained from the observations propagates 
both forward (to the right) and backward (Fig. 14b). How- 
ever, the forward propagation, which is in the direction of 
the group velocity (to the right), is faster than the back- 
ward propagation (to the left), and the maximum effect 
propagates at a speed that is close to the group velocity 
(section 7). 
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APPENDIX A 



Extracting an envelope 

See Zimin et al. (2003) for the original introduction of 
the method presented here. Let us say that we have a 
function 



f{n) = a{n) cos[ (l){n)], 0<n<N, 



(Al) 



where n is an integer, a(n) is an envelope, and cos[(j){n)] 
is an oscillating part with (f>{n) monotonically increasing. 
We want to extract a{n) given f{n). First, f{n) can be 
expressed as 



f{n) = ^a{n)e' 



(A2) 



Define F{k), A{k), E^{k), and E_{k) to be discrete Fourier 
transforms (DFTs) of f{n), a{n), e"''^"\ and e-*'^^"), re- 
spectively. Assuming each function is N-periodically ex- 
tended, we have the following: 



F{k) = ^DFT{a{n)e 



= —A{k)*E+{k) 



'} + i DFT{a(7i)e~ 
±^A{k)*E^{k) 



(A3) 
(A4) 



1 

— Y 

2N ^ 

1=0 



A{l)E+{k - I) 



N-l 

— Y 

2N ^ 

1=0 



A{l)E^{k~l), 
(A5) 



where the asterisk denotes circular convolution. The term 
E+(k) is centered around a positive frequency, and E-{k) 
is centered around a negative frequency. If we assume that 
a(n) changes much more slowly than 0(ri), then A{k) is well 
separated from E+{k) and E^{k). See Fig. 18, where the 
horizontal axis is the wavcnumber k and the vertical axis is 
the amplitude of the DFT values. The red, blue, and black 
plots are for DFTs of e^'^^"\ ^-^H^) ^ and a(ri), respectively. 
Therefore, A{k) * E+{k) will be on the positive frequency 
side and A(fc)*£'_ (fc) will be on the negative frequency side. 
So, taking inverse discrete Fourier transform of F{k) only 
with positive frequencies, DFT^^{i^(fc)}, is equivalent to 
taking inverse discrete Fourier transform of (2A^)^^A(fc) * 
E+{k). Therefore, we have the following: 



DFT:^^{F(A:)} = DFT" 



l-A{k)*E+{k) 



(A6) 
(A7) 
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Fig. 18. DFT of a(n), e*'^^"), and e-*'*^"). The terms ^(fc), 
E+{k), and E.{k) are DFTs of a(n), e"^("), and e-'^^"\ 
respectively, and f{n) = a{n) cos[(/)(n) ], Q <n < N. 



Thus, 

a(n)e^'^(") = 2DFT;^{F(fc)} (A8) 
Taking the absolute values of both sides, we obtain 



a(n) = 2|DFT;1{F(A:)}|. 



(A9) 



In the above equation, only several frequencies that have 
large DFT values can be used to approximate the enve- 
lope. This has the effect of smoothing out the envelope by 
sacrificing accuracy. As more and more frequency compo- 
nents are added, more and more wiggles are present in the 
envelope. 

APPENDIX B 



LETKF algorithm 

We present a summary of the LETKF algorithm of 
Hunt et al. (2007) used in our study with Lorenz model 
3. We assume a discrete model that has grid points. Each 
grid point has a variable. A local region associated with 
grid point n is defined as a region consisting of grid points 
i with n — L<i<n + L for a certain integer L. 

Input: a global background ensemble of TOg-dimensional 
model state vectors {x^'*'', i — 1,2, ... , k}, an Zg-dimensional 
vector y° of observations, a function H that maps the nig- 
dimensional model space to the ^^-dimensional observation 
space, and an Ig x Ig observation error covariance matrix 
Rg. Here the subscript g denotes global. 
Output: a global analysis ensemble of rrig-dimensional 

model state vectors {xg^'\ i = 1, 2, . . . , fc} 
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1. Apply H to each x^^'"* to form the global background 
observation ensemble {yg^'^} and average the latter vectors 
to get the Zg-dimensional column vector y^. Subtract this 

vector from each y^'-*'' to form the columns of the Ig x fc 
matrix Y^. 

2. Average the vectors {xg''-*} to get the rrig-dimcnsional 
column vector and subtract this vector from each Xg*-'-* 
to form the columns of the mg x k matrix X^. 

For each grid point, do the following steps 3-8. 

3. Select the rows of x^ and corresponding to the local 
region of m grid points associated with the given grid point, 
forming the m-dimensional vector x^ and the mx k matrix 
X^. Select the rows of y^ and corresponding to the I 
observations for the local region, forming the Z-dimensional 
vector y'' and the I x k matrix Y''. Select the corresponding 
rows of y°, forming the /-dimensional vector y°. Select the 
corresponding rows and columns of R^, forming the / x I 
matrix R. 

4. Compute the k x I matrix C (Y'')^R"\ 

5. Compute the kx k matrix P° = [{k - l)\/p + CY^]-\ 
where p > 1 is a covariance inflation factor. 

6. Compute the fc x fc matrix W = [(fc - 1)9°^]^/^, where 
the power 1/2 means the symmetric square root. 

7. Compute the fc-dimcnsional vector w° = P C(y° — y**) 
and add it to each column of W", forming a fc* x fc matrix 
whose columns are the analysis vectors {w"^'^}. 

8. Multiply X** by each w"(') and add x'' to get the local 
analysis ensemble members {x"*-*-*} for the given grid point. 

9. After performing steps 3-8 for each grid point, form 
the global analysis ensemble {xg''*'' } by combining 

from each local region by taking weighted averages, Xn[g = 
X]j=n-L ^n^/(j)/('^ ~ j)i where the subscript l{j) denotes 
the local region associated with the grid point j and /(•) is 
a weighting function defined on the integers from —L to L 
for a certain value L. 
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